BCR Work-Precision Diagrams
Samuel Isaacson and Chris Rackauckas
The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff chemical reaction network modeling the BCR signaling network from Barua et al.. We use ReactionNetworkImporters to load the BioNetGen model files as a Catalyst model, and then use ModelingToolkit to convert the Catalyst network model to ODEs.
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq, LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools gr() datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() tf = 100000.0 # generate ModelingToolkit ODEs @timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) rn = prnbng.rn @timeit to "Create ODESys" osys = convert(ODESystem, rn) tspan = (0.,tf) @timeit to "ODEProb No Jac" oprob = ODEProblem(osys, Float64[], tspan, Float64[]) @timeit to "ODEProb DenseJac" densejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true)
Parsing parameters...done
Creating parameters...done
Parsing species...done
Creating species...done
Creating species and parameters for evaluating expressions...done
Parsing and adding reactions...done
Parsing groups...done
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100000.0)
u0: 1122-element Vector{Float64}:
299717.8348854
47149.15480798
46979.01102231
290771.2428252
299980.7396749
300000.0
141.3151575495
0.1256496403614
0.4048783555301
140.8052338618
⋮
1.005585387399e-24
6.724953378237e-17
3.395560698281e-16
1.787990228838e-5
8.761844379939e-13
0.0002517949074779
0.0005539124513976
2.281251822741e-14
1.78232055967e-8
@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true, sparse=true) show(to)
──────────────────────────────────────────────────────────────────────────
──
Time Allocations
────────────────────── ─────────────────────
──
Tot / % measured: 406s / 99.0% 69.8GiB / 99.4%
Section ncalls time %tot avg alloc %tot a
vg
──────────────────────────────────────────────────────────────────────────
──
ODEProb DenseJac 1 280s 69.6% 280s 45.8GiB 66.0% 45.8G
iB
ODEProb SparseJac 1 62.4s 15.5% 62.4s 14.8GiB 21.3% 14.8G
iB
ODEProb No Jac 1 34.0s 8.45% 34.0s 4.78GiB 6.88% 4.78G
iB
Create ODESys 1 13.9s 3.47% 13.9s 2.61GiB 3.76% 2.61G
iB
Parse Network 1 11.7s 2.91% 11.7s 1.42GiB 2.05% 1.42G
iB
──────────────────────────────────────────────────────────────────────────
──
@show numspecies(rn) # Number of ODEs @show numreactions(rn) # Apprx. number of terms in the ODE @show length(parameters(rn)) # Number of Parameters
numspecies(rn) = 1122 numreactions(rn) = 24388 length(parameters(rn)) = 128 128
Time ODE derivative function compilation
As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.
u = ModelingToolkit.varmap_to_vars(nothing, species(rn); defaults=ModelingToolkit.defaults(rn)) du = copy(u) p = ModelingToolkit.varmap_to_vars(nothing, parameters(rn); defaults=ModelingToolkit.defaults(rn)) @timeit to "ODE rhs Eval1" oprob.f(du,u,p,0.) @timeit to "ODE rhs Eval2" oprob.f(du,u,p,0.) densejacprob.f(du,u,p,0.) sparsejacprob.f(du,u,p,0.)
We also time the ODE rhs function with BenchmarkTools as it is more accurate given how fast evaluating f is:
@btime oprob.f($du,$u,$p,0.)
38.769 μs (3 allocations: 512 bytes)
Now we time the Jacobian functions, including compilation time in the first evaluations
J = zeros(length(u),length(u)) @timeit to "DenseJac Eval1" densejacprob.f.jac(J,u,p,0.) @timeit to "DenseJac Eval2" densejacprob.f.jac(J,u,p,0.)
ERROR: syntax: expression too large
Js = similar(sparsejacprob.f.jac_prototype) @timeit to "SparseJac Eval1" sparsejacprob.f.jac(Js,u,p,0.) @timeit to "SparseJac Eval2" sparsejacprob.f.jac(Js,u,p,0.) show(to)
──────────────────────────────────────────────────────────────────────────
──
Time Allocations
────────────────────── ─────────────────────
──
Tot / % measured: 795s / 97.3% 86.6GiB / 99.3%
Section ncalls time %tot avg alloc %tot a
vg
──────────────────────────────────────────────────────────────────────────
──
SparseJac Eval1 1 303s 39.1% 303s 8.67GiB 10.1% 8.67G
iB
ODEProb DenseJac 1 280s 36.1% 280s 45.8GiB 53.3% 45.8G
iB
ODE rhs Eval1 1 67.1s 8.68% 67.1s 5.88GiB 6.84% 5.88G
iB
ODEProb SparseJac 1 62.4s 8.07% 62.4s 14.8GiB 17.2% 14.8G
iB
ODEProb No Jac 1 34.0s 4.39% 34.0s 4.78GiB 5.56% 4.78G
iB
Create ODESys 1 13.9s 1.80% 13.9s 2.61GiB 3.04% 2.61G
iB
Parse Network 1 11.7s 1.51% 11.7s 1.42GiB 1.65% 1.42G
iB
DenseJac Eval1 1 2.33s 0.30% 2.33s 2.01GiB 2.34% 2.01G
iB
SparseJac Eval2 1 82.7μs 0.00% 82.7μs 944B 0.00% 94
4B
ODE rhs Eval2 1 54.1μs 0.00% 54.1μs 688B 0.00% 68
8B
──────────────────────────────────────────────────────────────────────────
──
Picture of the solution
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5) plot(sol, legend=false, fmt=:png)
